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Abstract 

We test the applicability of density functional theory (DFT) to spectral perturbations taking 
CN . an example of a Cs atom surrounded by superfluid helium. The atomic DFT of helium is used 

to obtain the distribution of helium atoms around the impurity atom, and the electronic DFT 
is applied to the excitations of the atom, averaging over the ensemble of helium configurations. 
The shift and broadening of the D\ and D2 absorption lines are quite well reproduced by theory, 
suggesting that the DFT may be useful for describing spectral perturbations in more complex 
environments. 
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I. INTRODUCTION 



The time- dependent density-functional theory has proven to be a powerful tool in the 
description of optical absorption for molecules and clusters in free space (For a review, see 
Ref. [JJ. More recent citations can be found in Refs. J| [3|].). We would like to know whether 
the theory can be extended to calculate the line shifts and broadening when the absorber 
is embedded in a medium. Most applications have considered the electronic excitation from 
a single frozen nuclear configuration representing the ground state of an isolated molecular 
or an ideal crystal, in which case the predicted absorption lines below the photoionization 
threshold are sharp. Recently it was shown that the TDDFT also works well in describing 
the broadening of the transitions and the strength of symmetry-forbidden transitions due 
to zero-point vibrational motion, taking the example of the benzene spectrum J|. In the 
present work, we will calculate the effects of external perturbations on the optical absorption. 
We choose as a simple test case a Cs atom immersed in liquid helium at low temperature, 
because of the simplicity both of the electronic structure and of the external perturbation. 
In the long term we are interested in extending this kind of analysis more complex systems 
which would require the full power of the time-dependent density functional theory. 

Aside from our motivation from the perspective of applications of DFT, spectroscopic 
measurements of impurity atoms and molecules in superfluid helium have been attracting 
considerable interest in recent years || . The repulsive force between an impurity and helium 
atoms induces a "bubble" around the impurity. This leads to a weak perturbation of helium 
atoms on the spectra of impurities. The line shifts and spectral shapes induced by the helium 
perturbation provide information on the properties of the bubble in the quantum liquid as 
well as the excited states of the impurity. Since the perturbation is weak, this method 
also provides a unique opportunity to spectroscopic measurements of atomic clusters at low 
temperature || 0, ||. 

The perturbations of Cs lines have been studied experimentally as a function of helium 
density Q and we shall calculate this system. There are two s-to-p transitions, the D\ 
( s i/2 — * P1/2) an d D 2 (si/2 — ■> P3/2) lines, which are blue-shifted in helium and acquire 
widths. The shifts and widths of the two lines are different, and the D2 line has a skewed 
shape suggesting a double-peak structure. These features were first analyzed with a collective 



vibration model of the helium bubble [ 10| . The model reproduced average peak shifts, 
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but gave line widths less than a half of observed ones. A more sophisticated analysis has 
been made treating the liquid helium environment by the Path-Integral Monte-Carlo |IT 



This quantum simulation has succeeded to reproduce observed D 2 line profiles of the Cs 
spectra. However, this computation method is very costly in computer resources. We are 
thus motivated to develop an approach to treat helium perturbation that is both simple 
and yet has quantitative accuracy. We will show that a density functional theory (DFT) 
together with a statistical description of helium configurations meets our purpose. The 
helium density distribution around an embedded atom is calculated with the DFT, and the 
helium configurations are generated by a random sampling using the density distribution as 
the sampling weight function. 

Besides reporting the calculations on the absorption spectra of Cs in helium, we offer some 
simple qualitative interpretation of how the qualitative features of the spectrum reflect the 
properties of the helium bubble around the atom. 

II. FORMALISM 

A. Description of liquid helium around an impurity atom 

Among number of density functional methods for liquid helium, we adopt the Orsay- 
Paris functional of Ref. [|12[]. Although the Orsay- Paris functional is known to have some 
problems with dynamic properties of liquid helium Jl3| , it has a correct long-range behavior 



and reasonable short-range characteristics. Since we are interested in a density profile of 
liquid helium, it should be adequate for our purposes. 
The energy in the DFT is assumed to have the form 

E = j drH (r) } (1) 
vVpWf + \J rfrV(r)p(r')V LJ (|r - r'|) + °-p{v) (p r ) 1+ ^ . (2) 



where 



Here, m is the mass of a helium atom and p r is a coarse-grained density defined by 
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The VLj is a standard Lennard- Jones potential describing the He-He interaction screened at 
distances shorter than the distance h, 



46 











(]r-r'l) 





for |r — r'| > h, 

(4) 



VLj(h)(£p) , for \r-r'\<h. 



VLj(|r-r'|) 



The values of the parameters in Eqs. (@), (|), and (§ are c = 1.04554 x 10 7 KA 3 ( 1+7 ), 
7 = 2.8, e = 10.22 K, a = 2.556 A, and h = 2.377 A. This is the same density functional 



that was used in Ref. [T^ to study atomic impurities in liquid helium. In that work, the 
effect of the impurity was treated by including a potential interaction Vj-(r) between it and 
the helium atoms in the density functional, 

W(r)=Wo(r) + Fi(r)p(r), (5) 

The Vj(r) has important contributions from the repulsive between electrons and the helium 
atoms, as well as the Van der Waals-type polarization interaction. Since we need to treat 
interaction of impurity electrons with the helium atoms explicitly later on when we calculate 
the electronic excitation, we introduce it here as well for calculating the helium distribution. 
We approximate it as a contact interaction, i.e. of the form 

K-He(r e -r) = V q 6(t -t), (6) 

where r e and r are coordinates of the electron and helium atom, respectively. The strength 
Vo is determined from the electron-helium scattering length a as 

V = (7) 

m e 

where m e is the electron mass. Then Vi(r) is given by 

V x {v) = V oPc (r). (8) 

where p e ( r ) is the electron density of the impurity atom. We take a = 0.69 A, corresponding 
the observed low-energy electron- helium cross section a = 6.0 A 2 |TJ|]. We have also assumed 



that the ion core is heavy enough to be treated as a classical particle at the origin. Since 
Eq. (|6|) expresses the interaction between He atoms and an electron, the same interaction 
will be used to estimate the energy shift of valence electrons due to the helium perturbation. 
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This treatment of interatomic potential ignores long-range attraction due to the polarization 
effects. The influence of the polarization effect will be mentioned later. 

Utilizing the energy functional, E[p] = J drTC(r), we calculate the density profile of liquid 
helium, putting the impurity atom at the origin. Minimizing the grand potential at zero 
temperature, Q = E[p{r)] — pN, leads to a Hartree-type equation 

f 



■V + U(r) + V,(r) 



= l*Sp(?), (9) 



where 



U(t) = J dv'p(r')V LJ (\r-r'\) + ^(p r y +1 + £(1+ 7 )JL J ^dv' p{v>) {p r ,y . (fO) 

The equation is solved with the boundary condition that the density go to the bulk density 
Po at large r, which can be satisfied setting the chemical potential to 

/i=6 P o+(i+|)c P r i , (ff) 

where b = J drVLj(r) = — 8.888f x 10 2 KA 3 . The bulk density is related to the pressure P 
by 

Carrying out the solution of Eq. (H) we find the density profile shown in Fig. |l| The three 
curves give p(r) at equilibrium densities p of 0.0218, 0.0239, and 0.0253 A~ 3 , corresponding 
to P = 0, 10, and 20 atm, respectively. One can see a sharp rise in the density at r w 6 
A. This corresponds to the bubble radius. An oscillatory structure appears on the density 
profile, especially under high pressure. This feature is different from that of a bubble model 



adopted in Refs. 0, |10|- The maximum value of the density for P = 20 atm is about 0.0275 
A" 3 at r = 7.2 A. 

We use the p(r) computed above to generate configurations of helium atoms as follows. 
Take a large volume surrounding the alkali atom and denote it as V. This volume includes iV 
helium atoms in average, where iV is given by J v drp(r) = N. We randomly sample iV helium 
positions in V according to the density distribution p(r). This sampling procedure gives 
probability distribution without correlation among helium atoms. Denoting /(r) = p(r)/N, 
the probability distribution of iV atoms is given by 

N 

w nc (r)=Uf(r l ), (13) 
i=i 
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FIG. 1: Helium density profile around Cs atom. The coordinate r represents the distance from 
the Cs atom. 

where r stands for (r 1; • • • , rjv). We will also study the effect of helium- helium correlations 
by considering a probability distribution of the form 

N N 

w c (r)=l[g(r l )l[d(r l] -d). (14) 

i=l i<j 

Here d is the range of a short-range correlation. The distribution function g(r) is determined 
by the condition that the distribution of Eq. (|l~4l) gives a helium density p(r) , 

p(r) = N dr 2 ■ ■ -dr N w c {r,r 2 , ■ ■ ■ ,r N ). (15) 
Jv 

In practice, we employ an iterative procedure to find g(r) from this condition. 



B. Helium perturbation on the atomic spectra 

In the previous section, we described a density functional theory for calculating the 
density profile of liquid helium p(r), and taking account of the effect of an impurity atom 
at the origin. In this section, we discuss the calculation of the atomic spectrum, including 
the effects of the helium atoms. 

We begin with the theory of the isolated atom. Orbital wave functions are calculated using 
density functional theory with Dirac wave functions and kinetic energy operator. We need 
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accurate wave functions at large distances from the atom, which cannot be achieved with 
the traditional LDA functional due to the incorrect orbital eigenvalues and the incorrect 
asymptotic behavior of the potential. As is well known, these problems are diminished 
with the GGA density functionals. We employ the GGA functional of Ref. JIBf , which 
has a gradient correction to produce correct asymptotic behavior of the potential. The 
gradient correction includes an adjustable parameter j3\ we utilize this freedom to make the 
orbital energy coincide with the measured one. For the S1/2 orbital, the orbital energy is 
set equal to the ionization potential of Cs atom, 3.89 eV. For the py 2 and p 3 / 2 orbitals, 
the orbital energies are set equal to the ionization potential minus the excitation energies 
(about 1.43 eV). The quality of the wave function may be tested by examining the transition 
oscillator strength. For Di and D 2 transitions, the calculated oscillator strength assuming 
a pure single-electron transition is 1.034, in good agreement with measured value, 1.058. 
The calculated electron density distributions are shown in Fig. 0. In principle, there will be 
contributions to the transition from core electrons as well that can be taken into account with 
the TDDFT. Applying the TDDFT to the present case, the core contributions reduce the 
oscillator strength by some tens of percent, but do not significantly affect the asymptotic 
wave function of the valence electron. We therefore use the simpler single-electron wave 
functions below rather than wave functions from the TDDFT. 




0.001 



FIG. 2: Electron density distribution of the S]/ 2 (solid) and pi/ 2 orbitals of Cs. See text for details 
of the calculation. 
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For a given configuration of helium atoms, r = (ri, • • • , rjv), we evaluate shifts of orbital 
energies with first order perturbation theory. We use the same helium configuration for the 
ground state Si/ 2 and excited states pi/ 2 and P3/2, following the Frank-Condon principle. 
For S1/2 and pi/2 states, the energy shifts of the valence electron is calculated as 

A£«(r) = (^\Y;V c -KM^ k) ) = ^ W W (r*)| 2 ( 16 ) 

i i 

where k stands for orbital quantum numbers (£j) and either m state may be taken. For p 3 / 2 
states, the matrix elements depend on m and we have to diagonalize a 4 x 4 matrix to get 
the energy shifts. We then obtain two eigenenergies, each of which is doubly degenerate. 

Each helium configuration produces an energy shift and possible splitting but the tran- 
sitions remain sharp. The line broadening comes from the ensemble average over helium 
configurations. The line shape of the spectra for Di (si/2 —> P1/2) transition is given by 

S Dl (E) = [ drw{T)5 (E - (AE (p ^\t) - AE {s ^ 2 \ T ))) , (17) 
Jv 

where E is a shift from the energy position of the free atom. For the D2 (si/2 — > P3/2) 
transition, we need to add the two eigenmodes, 

S D2 (E) = j drw{r) {<5 (e - (a^ P3/2) (t) - A^/ 2 »(t))) 

+5 (E - (aE [ ^ /2 \t) - AEM( T ))Jj . (is) 

C. Origin of peak shifts, broadening, and splitting 

In this section, we discuss qualitatively argument how the blue shifts, the broadening 
of the lines, and the splitting of D 2 transitions occur. First consider a valence electron of 
the alkali atom interacting with a single helium atom at a position R with respect to the 
impurity atom. The energy shift of the valence electron is determined by the electron wave 
function at the position of helium atom, ^( r e = R)- 

For S1/2 and px/2 states, the energy shift Ae is calculated as 

Ae« = V oP i k \R), (19) 

as the first-order perturbation by K-He- Since Vo is positive, the energy shift, Eq. (0), is 
also positive. As may be seen from Fig. 0, the wave function of the pi/ 2 state is considerably 
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larger than the si/ 2 state outside the bubble. Thus, we expect a blue shift of Di excitation 
spectra, Ae (pl /2) - Ae (s V2) > 0. 

For p 3 / 2 states, the situation is slightly more complicated, because, in general, there are 
off-diagonal matrix elements among degenerate states with different m. However, all the off- 
diagonal elements vanish if we assume that the helium atom lies on the Z-axis, X = Y = 
and Z = R. We lose no generality in the case of a single helium atom. The helium atom at 
Z = R produces an energy shift Ae as 



\ , (/^) |r (,.:,,)> I 2Vbp? 3/2) , for |m| = l/2, 

Ae m = (ipm \V e -He\iPm ) = \ (20) 

0, for \m\ = 3/2. 



where is the angle-average electron density. There is no shift for \m\ = 3/2 states in the 
first-order perturbation by V c ^u c - As a result, the D 2 transition splits into two peaks; one 
has a blue shift, Ae^L - Ae (s V2) > o, and the other has a small red shift, -Ae (s V2) < o. 
In fact we can neglect — Ae^ 1 / 2 ^ in comparison to the other shifts due to the large distance 
of the perturbing helium atom. Also one can neglect the small energy difference between 
the two spin orbit partners. Then the energy shift of \m\ = 1/2 states is twice of that of 
P1/2 states. Note that the average shift of the D 2 components is the same as the D\ line. 
In summary, a helium atom at a distance R shifts the D\ transition by A = Vopf 1/2 \R). 
The energy of D 2 transition splits into two: one has no shift and the other is shifted by 2A. 
This explains qualitatively the overall blue shift of both lines and the skewed profile of the 
D 2 line. 

However, while the observed D 2 line shape can be analyzed as a sum of two components, 
both components are blue shifted roughly the same amount as the D\ line. This shows the 
shift is due to the simultaneous interaction with several helium atoms. For example, if two 
helium atoms are at (R, 0, 0) and at (0,0, R), again, a D 2 line splits into two. However, 
in this case, both have blue shifts (A and 3A). If three helium atoms are at (R, 0,0), 
(0, R, 0), and (0,0, R), then, D 2 line shows a blue shift of 3A but no splitting. Roughly 
speaking, the magnitude of line shift is determined by number of helium atoms contributing 
to the perturbation and the splitting is determined by anisotropy of the helium configuration 
which increases as the square root of the number of atoms, assuming that there are no 
correlations between atoms. The average energy shift remains the same for the D x and D 2 
lines, irrespective of the number of helium atoms causing the perturbation. 
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III. NUMERICAL RESULTS 



To calculate line shapes of the Cs D transition in a helium, we evaluated Eq. (2.18) 
and (2.19) by sampling 100 000 helium configurations, generated according to the density 
profiles in Fig. |1| The calculated energy shifts are added to the observed D lines of free 
Cs atom (A = 894.9 nm for D\ and 852.7 nm for D2). Then, the intensity is estimated 
by counting number of events in bins of wavelength AA = 0.1 nm. The obtained intensity 
spectra are shown in Fig. |3|. The D\ line can be well approximated by a single Gaussian, 
although there is a slightly larger tail at high-energy (low-wavelength) side. On the other 
hand, the D 2 line has a double-peaked structure. This feature agrees with the experimental 
observation |9|, [T(| . 



Fig. § shows the pressure dependence of the peak shift and broadening of the D\ excitation 
lines. The peak shift reproduces the observed pressure dependence, but comes out about 
20% lower than measured. The difference could be due to an incorrect asymptotic wave 
function in the Cs atom; only a 10% error in the wave function would be required to explain 
the difference. Or the calculated helium bubble might be too large. Here, decreasing the 
size of the bubble by 0.3 A out of 6 A would be sufficient to produce the measured peak 
shift. We shall return to this later. The line broadening comes out better than would be 
expected, given the quality of agreement for the peak shift. The agreement here shows that 
the fluctuations of the helium distribution are well described by the model adopted. 

A qualitative measure of the fluctuations can be constructed by defining an effective 
number of helium atoms that contribute to the perturbation. Calling the shift from an 
individual helium atom A(i), the total shift is 

N 

AE = Y J ^) (21) 
i=i 

and the effective number M e g is defined as 

M _/ (E,A(i)) 2 \ 

where (• • • ) indicates the ensemble average. The calculation leads to M e fj « 8 for syz state 
at P = 0, increasing to 11-12 at P = 25 atm. The M e g can be compared with number of 
helium atoms in the first shell of density profile (Fig. [l]). The average numbers of helium 
atoms in a region of r < 8.5 A are 27 at P = and 45 at P = 25 atm. M e g turns out 
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FIG. 3: (a) Cs D\ excitation spectra at different helium pressure; P = 0, 10, and 20 atm (b) The 
same as (a) but for D 2 excitation spectra. Experimental data from Ref. || are plotted as filled 
squares. 

to be much smaller than the number in the first shell. Therefore, we may say that the 
perturbation on the valence electron is dominated by a small number of helium atoms in the 
inner surface of a bubble. This fact indicates the importance of treating the perturbation 
from individual helium atoms in describing fluctuation effect. 

We can also understand the order of magnitude of the fluctuation effects using M e ff. 



Assuming independent helium atoms, the fluctuations are proportional to l/y/M c s. Thus 
we would expect that widths of the lines and the splitting of the components of the D 2 would 
be proportional to the average shift times that quantity. In fact the D\ width is about 1/2 
of its average shift in the zero pressure data, to be compared with 1 jyj M eff ps 1/3. Ref. 
H analyzes the D 2 line shape as a sum two Gaussian peaks. For the measurements at low 
pressure, the splitting of the components relative to the average shift is just 1/3, agreeing 
with our very crude argument. 

In order to investigate effects of the polarization potential that we so neglected, we 



consider also a helium-Cs interaction that includes the van der Waals terms 17 . Once the 
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FIG. 4: Pressure shift (solid line) and broadening (dashed line) of the D\ excitation spectra of 
Cs as a function of helium pressure. Open squares represent broadening calculated with He-He 
correlations at P = 0, 10, and 20 atm. The experimental data are plotted as filled squares. 

density profile p(r) is determined, we use the interaction between a valence electron and 
helium, Eq. @, to calculate the atomic spectra. We find that potential of Ref. |l7j gives 
a reduced radius for the helium bubble at zero pressure, by about 0.3 A. As a result, the 
blue shift increases by about 35 cm -1 at P = 0, which fits the experimental data very well. 
The line width is also slightly increased. This enhancement of blue shifts is decreasing as 
the pressure is increasing. The calculation shows almost no additional shift at P=20 atm. 

Next, let us discuss effects of correlations among helium atoms. The He- He correlation 
should influence the line width because it removes some part of fluctuations of the He 
configuration. Roughly speaking, the radial fluctuation controls the line width and the 
angle fluctuation determines the skewness of D 2 line. To test this, we sample the helium 
configurations using the probability distribution of Eq. (0), taking d = 2.377 A. We 
construct the distribution function g(r) in Eq. (|T4]), so as to reproduce the density profile 
p(r) determined by the DFT, Eq. @. The results are displayed in Fig. [5] and also in Fig. |] 
with open squares. The correlation effect does not change the average peak positions at 
all. However, it reduces the line widths, especially when the liquid helium is under high 
pressure. At P = 20 atm, the FWHM of D\ line is calculated to be 150 cm -1 , which well 
agrees with the experiment (140 cm -1 ). The skewness of D 2 line also becomes smaller. 
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FIG. 5: (a) Cs D\ excitation spectra and (b) D2 excitation spectra as in Fig. 2, with inclusion of 
He-He correlations. See text for explanation. 

IV. SUMMARY 

We have developed a simple model to describe atomic spectra of impurities embedded 
in the superfluid helium. Our description employs a density functional theory for the he- 
lium distribution, and treats helium configurations statistically. The model is applied to 
the spectra of Cs atom embedded in superfluid helium. Various features in the spectra, 
which include line shifts, broadening, and skewness, are nicely reproduced in our calculation 
without any adjustable parameters. Thus we are confident that our model includes basic 
physical elements of the helium perturbation correctly. Since the model is simple enough 
to apply to more complex systems such as molecules and clusters, we wish to analyze these 
systems in future. 
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